Part1

library(spatstat)
library(here)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
library(tmap)
library(sf)
library(geojson)
library(geojsonio)
library(tmaptools)
library(stringr)
library(tidyverse)
library(spdep)
library(raster)
library(fpc)
library(dbscan)
library(OpenStreetMap)
library(png)
library(rmarkdown)
#reading wards
LondonWards <- st_read(here::here("data", "London-wards-2018_ESRI", "London_Ward.shp"))

## Reading layer `London_Ward' from data source `D:\OneDrive\UCL\Work\GIS\GIS-coursework\Data\Cutural infrastructures\data\London-wards-2018_ESRI\London_Ward.shp' using driver `ESRI Shapefile'
## Simple feature collection with 657 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS:  OSGB 1936 / British National Grid

LondonWardsMerged <- st_read(here::here("data", "London-wards-2018_ESRI",
                                        "London_Ward_CityMerged.shp"))%>%
  st_transform(.,27700)

## Reading layer `London_Ward_CityMerged' from data source `D:\OneDrive\UCL\Work\GIS\GIS-coursework\Data\Cutural infrastructures\data\London-wards-2018_ESRI\London_Ward_CityMerged.shp' using driver `ESRI Shapefile'
## Simple feature collection with 633 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS:  OSGB 1936 / British National Grid

#reading pubs
Pubs <- read_csv("data/Pubs.csv",na = c("NA", "n/a"))
Pubs_spatial <- sf::st_as_sf(Pubs, coords = c("longitude", "latitude"), crs = 4326) %>%
  st_transform(.,27700)
  
qtm(LondonWardsMerged)

qtm(Pubs_spatial)

#change the ward name to code
#library(plyr) 
#df <- data.frame(foo=norm(1000)) 
#df <- rename(df,c('foo'='samples'))
names(Pubs_spatial)[names(Pubs_spatial) == 'ward_2018_code'] <- 'GSS_CODE'
names(Pubs_spatial)[names(Pubs_spatial) == 'ward_2018_name'] <- 'NAME'
#join data

#LondonWardsMerged_Pub <- LondonWardsMerged
Pubs_spatial <- Pubs_spatial %>%
  add_count(GSS_CODE, name="Pubs_in_ward")
#don't know how to join the don't Let's do it
Pubs_sub <- Pubs_spatial[LondonWardsMerged,]
tmap_mode("view")
tm_shape(LondonWardsMerged) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(Pubs_sub) +
  tm_dots(col = "blue")

#mapping the density of points
library(sf)
LondonWardsMerged %>% 
  mutate(ward_area = st_area(geometry)) %>% 
  head()

## Simple feature collection with 6 features and 7 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 516362.6 ymin: 159907.4 xmax: 522194.7 ymax: 172367
## projected CRS:  OSGB 1936 / British National Grid
##                     NAME  GSS_CODE             DISTRICT LAGSSCODE HECTARES
## 1      Chessington South E05000405 Kingston upon Thames E09000021  755.173
## 2 Tolworth and Hook Rise E05000414 Kingston upon Thames E09000021  259.464
## 3             Berrylands E05000401 Kingston upon Thames E09000021  145.390
## 4              Alexandra E05000400 Kingston upon Thames E09000021  268.506
## 5               Beverley E05000402 Kingston upon Thames E09000021  187.821
## 6            Coombe Hill E05000406 Kingston upon Thames E09000021  442.170
##   NONLD_AREA                       geometry     ward_area
## 1          0 POLYGON ((516401.6 160201.8... 7551784 [m^2]
## 2          0 POLYGON ((519553 164295.6, ... 2594673 [m^2]
## 3          0 POLYGON ((518107.5 167303.4... 1453889 [m^2]
## 4          0 POLYGON ((520336.7 165105.5... 2685017 [m^2]
## 5          0 POLYGON ((521201.2 169275.5... 1878205 [m^2]
## 6          0 POLYGON ((521296.5 172258.3... 4421680 [m^2]

LondonWardsMerged %>% 
  mutate(ward_area = st_area(geometry)) %>% 
  st_join(Pubs_spatial) %>%
  head()

## Simple feature collection with 6 features and 21 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 516362.6 ymin: 159907.4 xmax: 520336.7 ymax: 168160.9
## projected CRS:  OSGB 1936 / British National Grid
##                     NAME.x GSS_CODE.x             DISTRICT LAGSSCODE HECTARES
## 1        Chessington South  E05000405 Kingston upon Thames E09000021  755.173
## 1.1      Chessington South  E05000405 Kingston upon Thames E09000021  755.173
## 1.2      Chessington South  E05000405 Kingston upon Thames E09000021  755.173
## 1.3      Chessington South  E05000405 Kingston upon Thames E09000021  755.173
## 2   Tolworth and Hook Rise  E05000414 Kingston upon Thames E09000021  259.464
## 3               Berrylands  E05000401 Kingston upon Thames E09000021  145.390
##     NONLD_AREA     ward_area           name             address1
## 1            0 7551784 [m^2]           Star     Leatherhead Road
## 1.1          0 7551784 [m^2] William Bourne       273, Moor Lane
## 1.2          0 7551784 [m^2]  Monkey Puzzle     Leatherhead Road
## 1.3          0 7551784 [m^2]      Shy Horse 423 Leatherhead Road
## 2            0 2594673 [m^2]           <NA>                 <NA>
## 3            0 1453889 [m^2]     Berrylands  108, Chiltern Drive
##               address2 address3         borough_name website
## 1   Leatherhead Common          Kingston upon Thames      NA
## 1.1        Chessington          Kingston upon Thames      NA
## 1.2        Chessington          Kingston upon Thames      NA
## 1.3     Malden Rushett          Kingston upon Thames      NA
## 2                 <NA>     <NA>                 <NA>      NA
## 3             Surbiton          Kingston upon Thames      NA
##     os_addressbase_uprn borough_code            NAME.y GSS_CODE.y easting
## 1                    NA    E09000021 Chessington South  E05000405  516633
## 1.1                  NA    E09000021 Chessington South  E05000405  519187
## 1.2                  NA    E09000021 Chessington South  E05000405  517459
## 1.3                  NA    E09000021 Chessington South  E05000405  517170
## 2                    NA         <NA>              <NA>       <NA>      NA
## 3                    NA    E09000021        Berrylands  E05000401  519726
##     northing    runtime Pubs_in_ward                       geometry
## 1   160000.9 01/03/2021            4 POLYGON ((516401.6 160201.8...
## 1.1 163995.9 01/03/2021            4 POLYGON ((516401.6 160201.8...
## 1.2 162405.9 01/03/2021            4 POLYGON ((516401.6 160201.8...
## 1.3 161293.9 01/03/2021            4 POLYGON ((516401.6 160201.8...
## 2         NA       <NA>           NA POLYGON ((519553 164295.6, ...
## 3   168001.9 01/03/2021            5 POLYGON ((518107.5 167303.4...

Pubs_density <- LondonWardsMerged %>% 
  mutate(ward_area = st_area(geometry)) %>%
  st_join(Pubs_spatial) %>%
  group_by(GSS_CODE.x) %>% 
  summarize(n_Pubs = n(),
            ward_area = unique(ward_area),
            pubsdensity = n_Pubs/ward_area * 1e6)
plot(Pubs_density["pubsdensity"])

png('points density.png')
dev.off()

## png 
##   2

#plot density map
breaks1 <- c(0,0.10,0.65,1.5,3,6.5,15,25,35,50,Inf)
tmap_mode("plot")
tm_shape(Pubs_density) +
  tm_polygons("pubsdensity", 
              style="fixed",
              palette = "OrRd",
              breaks=breaks1,
              title="London \nPubs\nDensity \nper sqKm") +
  tm_scale_bar(position=c("right", "bottom"))

#now for the Moran's I and other statistics
coordsW <- Pubs_density %>%
  st_centroid()%>%
  st_geometry()
plot(coordsW,axes=TRUE)

LWard_nb <- Pubs_density %>%
  poly2nb(., queen=T)
#plot them
plot(LWard_nb, st_geometry(coordsW), col="red")
#add a map underneath, not very clear about the matrix weight part
plot(Pubs_density$geometry, add=T)

Lward.lw <- LWard_nb %>%
  nb2listw(., style="C")

head(Lward.lw$neighbours)

## [[1]]
## [1]   6   7  10 380 386 597
## 
## [[2]]
## [1]  5  8  9 11 12 13 16
## 
## [[3]]
## [1]  10  11  12  15 594 598
## 
## [[4]]
## [1]  17 217 227 584 587 595
## 
## [[5]]
## [1]   2   9  16 217 219 226
## 
## [[6]]
## [1]  1  7  8 10 11 14

I_LWard_Global_Density <- Pubs_density %>%
  pull(pubsdensity) %>%
  as.vector()%>%
  moran.test(., Lward.lw)

I_LWard_Global_Density

## 
##  Moran I test under randomisation
## 
## data:  .  
## weights: Lward.lw    
## 
## Moran I statistic standard deviate = 32.642, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.7353283210     -0.0015822785      0.0005096696

#use the localmoran function to generate I for each ward in the city
I_LWard_Local_count <- Pubs_density %>%
  pull(n_Pubs) %>%
  as.vector()%>%
  localmoran(., Lward.lw)%>%
  as_tibble()

I_LWard_Local_Density <- Pubs_density %>%
  pull(pubsdensity) %>%
  as.vector()%>%
  localmoran(., Lward.lw)%>%
  as_tibble()
#what does the output (the localMoran object) look like?
slice_head(I_LWard_Local_Density, n=5)

## # A tibble: 5 x 5
##   Ii          E.Ii         Var.Ii     Z.Ii        `Pr(z > 0)`
##   <localmrn>  <localmrn>   <localmrn> <localmrn>  <localmrn> 
## 1 -0.01181041 -0.001611124 0.1642175  -0.02516866 0.5100398  
## 2  0.29599684 -0.001879645 0.1912947   0.68105899 0.2479171  
## 3  0.25008892 -0.001611124 0.1642175   0.62111717 0.2672613  
## 4  0.20663155 -0.001611124 0.1642175   0.51387796 0.3036687  
## 5  0.27205882 -0.001611124 0.1642175   0.67533203 0.2497324

#copy some of the columns to the LondonWards spatialPolygonsDataframe
Pubs_density <- Pubs_density %>%
  mutate(pubs_count_I = as.numeric(I_LWard_Local_count$Ii))%>%
  mutate(pubs_count_Iz =as.numeric(I_LWard_Local_count$Z.Ii))%>%
  mutate(density_I =as.numeric(I_LWard_Local_Density$Ii))%>%
  mutate(density_Iz =as.numeric(I_LWard_Local_Density$Z.Ii))

#now plot moran's I
tmap_mode("plot") 
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
MoranColours<- rev(brewer.pal(8, "RdGy"))
tm_shape(Pubs_density) +
  tm_polygons("pubs_count_Iz",
              style="fixed",
              breaks=breaks1,
              palette=MoranColours,
              midpoint=NA,
              title="Local Moran's I, Pubs in London")

# now for the Getis Ord G
Gi_LWard_Local_Density <- Pubs_density %>%
  pull(pubsdensity) %>%
  as.vector()%>%
  localG(., Lward.lw)
head(Gi_LWard_Local_Density)

## [1] -1.1788424 -1.3311527 -1.2280027 -1.2499370 -1.2688086 -0.9549045

Pubs_density <- Pubs_density %>%
  mutate(density_G = as.numeric(Gi_LWard_Local_Density))
GIColours<- rev(brewer.pal(8, "RdBu"))

#now plot 
tm_shape(Pubs_density) +
  tm_polygons("density_G",
              style="fixed",
              breaks=breaks1,
              palette=GIColours,
              midpoint=NA,
              title="Gi*, Pubs in London")